Galaxy Liver Mets Kataoka et al_07082024 Clinical analysis

library(swimplot) library(coxphf) library(grid) library(gtable) library(readr) library(mosaic) library(dplyr) library(survival) library(survminer) library(ggplot2) library(scales) library(ggthemes) library(tidyverse) library(gtsummary) library(flextable) library(parameters) library(car) library(grid) library(ComplexHeatmap) library(readxl) library(janitor) library(rms) library(DT)

#Demographics Table

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]

circ_data_subset <- circ_data %>%
  select(
    Age,
    Gender,
    ECOG,
    PrimSite,
    Mets.Type,
    Hepatectomy,
    NumLiverMetsGroup,
    SizeLiverMetsmmGroup,
    pT,
    pN,
    ACT,
    Oxaliplatin.History,
    Postop.Complication,
    BRAF.V600E,
    RAS,
    MSI,
    RFS.Event,
    OS.months) %>%
  mutate(
    Age = as.numeric(Age),
    Gender = factor(Gender, levels = c("Male", "Female")),
    ECOG = factor(ECOG, levels = c(0, 1)),
    PrimSite = factor(PrimSite, levels = c("Right-sided colon", "Left-sided colon")),
    Mets.Type = factor(Mets.Type, levels = c("Synchronous", "Metachronous")),
    Hepatectomy = factor(Hepatectomy, levels = c("Minor", "Major")),
    NumLiverMetsGroup = factor(NumLiverMetsGroup, levels = c("1", "≥2")),
    SizeLiverMetsmmGroup = factor(SizeLiverMetsmmGroup, levels = c("<50", "≥50")),
    pT = factor(pT, levels = c("T1-T2", "T3-T4")),
    pN = factor(pN, levels = c("N0", "N1-N2")),
    ACT = factor(ACT, levels = c("TRUE", "FALSE"), labels = c("Adjuvant Chemotherapy", "Observation")),
    Oxaliplatin.History = factor(Oxaliplatin.History, levels = c("No", "Yes")),
    Postop.Complication = factor(Postop.Complication, levels = c("No", "Yes")),
    BRAF.V600E = factor(BRAF.V600E, levels = c("WT", "MUT"), labels = c("BRAF wt", "BRAF V600E")),
    RAS = factor(RAS, levels = c("WT", "MUT"), labels = c("RAS wt", "RAS mut")),
    MSI = factor(MSI, levels = c("MSS", "MSI-High")),
    RFS.Event = factor(RFS.Event, levels = c("TRUE", "FALSE"), labels = c("Recurrence", "No Recurrence")),
    OS.months = as.numeric(OS.months))
table1 <- circ_data_subset %>%
  tbl_summary(
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)")) %>%
  bold_labels()
table1

Characteristic	N = 1901
Age	68 (34 - 85)
Gender
Male	118 (62%)
Female	72 (38%)
ECOG
0	186 (98%)
1	4 (2.1%)
PrimSite
Right-sided colon	50 (26%)
Left-sided colon	140 (74%)
Mets.Type
Synchronous	64 (34%)
Metachronous	126 (66%)
Hepatectomy
Minor	159 (84%)
Major	31 (16%)
NumLiverMetsGroup
1	121 (64%)
≥2	69 (36%)
SizeLiverMetsmmGroup
<50	180 (95%)
≥50	10 (5.3%)
pT
T1-T2	23 (12%)
T3-T4	167 (88%)
pN
N0	77 (41%)
N1-N2	113 (59%)
ACT
Adjuvant Chemotherapy	48 (25%)
Observation	142 (75%)
Oxaliplatin.History	50 (26%)
Postop.Complication	23 (12%)
BRAF.V600E
BRAF wt	188 (99%)
BRAF V600E	2 (1.1%)
RAS
RAS wt	106 (56%)
RAS mut	84 (44%)
MSI
MSS	188 (99%)
MSI-High	2 (1.1%)
RFS.Event
Recurrence	95 (50%)
No Recurrence	95 (50%)
OS.months	24 (1 - 48)
1 Median (Range); n (%)

fit1 <- as_flex_table(
  table1,
  include = everything(),
  return_calls = FALSE,
  strip_md_bold = TRUE)

Warning: The `strip_md_bold` argument of `as_flex_table()` is deprecated as of gtsummary 1.6.0.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

fit1

Characteristic	N = 1901
Age	68 (34 - 85)
Gender
Male	118 (62%)
Female	72 (38%)
ECOG
0	186 (98%)
1	4 (2.1%)
PrimSite
Right-sided colon	50 (26%)
Left-sided colon	140 (74%)
Mets.Type
Synchronous	64 (34%)
Metachronous	126 (66%)
Hepatectomy
Minor	159 (84%)
Major	31 (16%)
NumLiverMetsGroup
1	121 (64%)
≥2	69 (36%)
SizeLiverMetsmmGroup
<50	180 (95%)
≥50	10 (5.3%)
pT
T1-T2	23 (12%)
T3-T4	167 (88%)
pN
N0	77 (41%)
N1-N2	113 (59%)
ACT
Adjuvant Chemotherapy	48 (25%)
Observation	142 (75%)
Oxaliplatin.History	50 (26%)
Postop.Complication	23 (12%)
BRAF.V600E
BRAF wt	188 (99%)
BRAF V600E	2 (1.1%)
RAS
RAS wt	106 (56%)
RAS mut	84 (44%)
MSI
MSS	188 (99%)
MSI-High	2 (1.1%)
RFS.Event
Recurrence	95 (50%)
No Recurrence	95 (50%)
OS.months	24 (1 - 48)
1Median (Range); n (%)

save_as_docx(fit1, path= "~/Downloads/table1.docx")

#Demographics Table by MRD ctDNA Status

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]

circ_data_subset1 <- circ_data %>%
  select(
    Age,
    Gender,
    ECOG,
    PrimSite,
    Mets.Type,
    Hepatectomy,
    NumLiverMetsGroup,
    SizeLiverMetsmmGroup,
    pT,
    pN,
    ACT,
    Oxaliplatin.History,
    Postop.Complication,
    BRAF.V600E,
    RAS,
    MSI) %>%
  mutate(
    Age = as.numeric(Age),
    Gender = factor(Gender, levels = c("Male", "Female")),
    ECOG = factor(ECOG, levels = c(0, 1)),
    PrimSite = factor(PrimSite, levels = c("Right-sided colon", "Left-sided colon")),
    Mets.Type = factor(Mets.Type, levels = c("Synchronous", "Metachronous")),
    Hepatectomy = factor(Hepatectomy, levels = c("Minor", "Major")),
    NumLiverMetsGroup = factor(NumLiverMetsGroup, levels = c("1", "≥2")),
    SizeLiverMetsmmGroup = factor(SizeLiverMetsmmGroup, levels = c("<50", "≥50")),
    pT = factor(pT, levels = c("T1-T2", "T3-T4")),
    pN = factor(pN, levels = c("N0", "N1-N2")),
    ACT = factor(ACT, levels = c("TRUE", "FALSE"), labels = c("Adjuvant Chemotherapy", "Observation")),
    Oxaliplatin.History = factor(Oxaliplatin.History, levels = c("No", "Yes")),
    Postop.Complication = factor(Postop.Complication, levels = c("No", "Yes")),
    BRAF.V600E = factor(BRAF.V600E, levels = c("WT", "MUT"), labels = c("BRAF wt", "BRAF V600E")),
    RAS = factor(RAS, levels = c("WT", "MUT"), labels = c("RAS wt", "RAS mut")),
    MSI = factor(MSI, levels = c("MSS", "MSI-High")))

circ_data1 <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data1 <- circ_data1[circ_data1$LiverMets=="TRUE",]

circ_data_subset2 <- circ_data %>%
  select(
    Age,
    Gender,
    ECOG,
    PrimSite,
    Mets.Type,
    Hepatectomy,
    NumLiverMetsGroup,
    SizeLiverMetsmmGroup,
    pT,
    pN,
    ACT,
    Oxaliplatin.History,
    Postop.Complication,
    BRAF.V600E,
    RAS,
    MSI,
    ctDNA.MRD) %>%
  mutate(
    Age = as.numeric(Age),
    Gender = factor(Gender, levels = c("Male", "Female")),
    ECOG = factor(ECOG, levels = c(0, 1)),
    PrimSite = factor(PrimSite, levels = c("Right-sided colon", "Left-sided colon")),
    Mets.Type = factor(Mets.Type, levels = c("Synchronous", "Metachronous")),
    Hepatectomy = factor(Hepatectomy, levels = c("Minor", "Major")),
    NumLiverMetsGroup = factor(NumLiverMetsGroup, levels = c("1", "≥2")),
    SizeLiverMetsmmGroup = factor(SizeLiverMetsmmGroup, levels = c("<50", "≥50")),
    pT = factor(pT, levels = c("T1-T2", "T3-T4")),
    pN = factor(pN, levels = c("N0", "N1-N2")),
    ACT = factor(ACT, levels = c("TRUE", "FALSE"), labels = c("Adjuvant Chemotherapy", "Observation")),
    Oxaliplatin.History = factor(Oxaliplatin.History, levels = c("No", "Yes")),
    Postop.Complication = factor(Postop.Complication, levels = c("No", "Yes")),
    BRAF.V600E = factor(BRAF.V600E, levels = c("WT", "MUT"), labels = c("BRAF wt", "BRAF V600E")),
    RAS = factor(RAS, levels = c("WT", "MUT"), labels = c("RAS wt", "RAS mut")),
    MSI = factor(MSI, levels = c("MSS", "MSI-High")),
    ctDNA.MRD = factor(ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE")))
Overall <- circ_data_subset1 %>%
  tbl_summary(
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)")) %>%
  bold_labels()
Overall

Characteristic	N = 1901
Age	68 (34 - 85)
Gender
Male	118 (62%)
Female	72 (38%)
ECOG
0	186 (98%)
1	4 (2.1%)
PrimSite
Right-sided colon	50 (26%)
Left-sided colon	140 (74%)
Mets.Type
Synchronous	64 (34%)
Metachronous	126 (66%)
Hepatectomy
Minor	159 (84%)
Major	31 (16%)
NumLiverMetsGroup
1	121 (64%)
≥2	69 (36%)
SizeLiverMetsmmGroup
<50	180 (95%)
≥50	10 (5.3%)
pT
T1-T2	23 (12%)
T3-T4	167 (88%)
pN
N0	77 (41%)
N1-N2	113 (59%)
ACT
Adjuvant Chemotherapy	48 (25%)
Observation	142 (75%)
Oxaliplatin.History	50 (26%)
Postop.Complication	23 (12%)
BRAF.V600E
BRAF wt	188 (99%)
BRAF V600E	2 (1.1%)
RAS
RAS wt	106 (56%)
RAS mut	84 (44%)
MSI
MSS	188 (99%)
MSI-High	2 (1.1%)
1 Median (Range); n (%)

ByctDNA_MRD <- circ_data_subset2 %>%
  tbl_summary(
    by = ctDNA.MRD, # add this line to subgroup by ctDNA.MRD
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)")) %>%
  add_p() %>%
  bold_labels()
ByctDNA_MRD

Characteristic	NEGATIVE, N = 1291	POSITIVE, N = 611	p-value2
Age	69 (37 - 85)	68 (34 - 84)	0.4
Gender			0.4
Male	83 (64%)	35 (57%)
Female	46 (36%)	26 (43%)
ECOG			0.6
0	127 (98%)	59 (97%)
1	2 (1.6%)	2 (3.3%)
PrimSite			0.033
Right-sided colon	40 (31%)	10 (16%)
Left-sided colon	89 (69%)	51 (84%)
Mets.Type			<0.001
Synchronous	33 (26%)	31 (51%)
Metachronous	96 (74%)	30 (49%)
Hepatectomy			0.089
Minor	112 (87%)	47 (77%)
Major	17 (13%)	14 (23%)
NumLiverMetsGroup			<0.001
1	93 (72%)	28 (46%)
≥2	36 (28%)	33 (54%)
SizeLiverMetsmmGroup			0.7
<50	123 (95%)	57 (93%)
≥50	6 (4.7%)	4 (6.6%)
pT			0.9
T1-T2	16 (12%)	7 (11%)
T3-T4	113 (88%)	54 (89%)
pN			0.015
N0	60 (47%)	17 (28%)
N1-N2	69 (53%)	44 (72%)
ACT			0.9
Adjuvant Chemotherapy	33 (26%)	15 (25%)
Observation	96 (74%)	46 (75%)
Oxaliplatin.History	39 (30%)	11 (18%)	0.075
Postop.Complication	10 (7.8%)	13 (21%)	0.007
BRAF.V600E			>0.9
BRAF wt	127 (98%)	61 (100%)
BRAF V600E	2 (1.6%)	0 (0%)
RAS			0.5
RAS wt	74 (57%)	32 (52%)
RAS mut	55 (43%)	29 (48%)
MSI			0.5
MSS	128 (99%)	60 (98%)
MSI-High	1 (0.8%)	1 (1.6%)
1 Median (Range); n (%)
2 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test

merged_table <- tbl_merge(tbls=list(Overall, ByctDNA_MRD))
merged_table

Characteristic	Table 1	Table 2
	N = 1901	NEGATIVE, N = 1291	POSITIVE, N = 611	p-value2
Age	68 (34 - 85)	69 (37 - 85)	68 (34 - 84)	0.4
Gender				0.4
Male	118 (62%)	83 (64%)	35 (57%)
Female	72 (38%)	46 (36%)	26 (43%)
ECOG				0.6
0	186 (98%)	127 (98%)	59 (97%)
1	4 (2.1%)	2 (1.6%)	2 (3.3%)
PrimSite				0.033
Right-sided colon	50 (26%)	40 (31%)	10 (16%)
Left-sided colon	140 (74%)	89 (69%)	51 (84%)
Mets.Type				<0.001
Synchronous	64 (34%)	33 (26%)	31 (51%)
Metachronous	126 (66%)	96 (74%)	30 (49%)
Hepatectomy				0.089
Minor	159 (84%)	112 (87%)	47 (77%)
Major	31 (16%)	17 (13%)	14 (23%)
NumLiverMetsGroup				<0.001
1	121 (64%)	93 (72%)	28 (46%)
≥2	69 (36%)	36 (28%)	33 (54%)
SizeLiverMetsmmGroup				0.7
<50	180 (95%)	123 (95%)	57 (93%)
≥50	10 (5.3%)	6 (4.7%)	4 (6.6%)
pT				0.9
T1-T2	23 (12%)	16 (12%)	7 (11%)
T3-T4	167 (88%)	113 (88%)	54 (89%)
pN				0.015
N0	77 (41%)	60 (47%)	17 (28%)
N1-N2	113 (59%)	69 (53%)	44 (72%)
ACT				0.9
Adjuvant Chemotherapy	48 (25%)	33 (26%)	15 (25%)
Observation	142 (75%)	96 (74%)	46 (75%)
Oxaliplatin.History	50 (26%)	39 (30%)	11 (18%)	0.075
Postop.Complication	23 (12%)	10 (7.8%)	13 (21%)	0.007
BRAF.V600E				>0.9
BRAF wt	188 (99%)	127 (98%)	61 (100%)
BRAF V600E	2 (1.1%)	2 (1.6%)	0 (0%)
RAS				0.5
RAS wt	106 (56%)	74 (57%)	32 (52%)
RAS mut	84 (44%)	55 (43%)	29 (48%)
MSI				0.5
MSS	188 (99%)	128 (99%)	60 (98%)
MSI-High	2 (1.1%)	1 (0.8%)	1 (1.6%)
1 Median (Range); n (%)
2 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test

fit1 <- as_flex_table(
  merged_table,
  include = everything(),
  return_calls = FALSE,
  strip_md_bold = TRUE)
fit1

	Table 1	Table 2
Characteristic	N = 1901	NEGATIVE, N = 1291	POSITIVE, N = 611	p-value2
Age	68 (34 - 85)	69 (37 - 85)	68 (34 - 84)	0.4
Gender				0.4
Male	118 (62%)	83 (64%)	35 (57%)
Female	72 (38%)	46 (36%)	26 (43%)
ECOG				0.6
0	186 (98%)	127 (98%)	59 (97%)
1	4 (2.1%)	2 (1.6%)	2 (3.3%)
PrimSite				0.033
Right-sided colon	50 (26%)	40 (31%)	10 (16%)
Left-sided colon	140 (74%)	89 (69%)	51 (84%)
Mets.Type				<0.001
Synchronous	64 (34%)	33 (26%)	31 (51%)
Metachronous	126 (66%)	96 (74%)	30 (49%)
Hepatectomy				0.089
Minor	159 (84%)	112 (87%)	47 (77%)
Major	31 (16%)	17 (13%)	14 (23%)
NumLiverMetsGroup				<0.001
1	121 (64%)	93 (72%)	28 (46%)
≥2	69 (36%)	36 (28%)	33 (54%)
SizeLiverMetsmmGroup				0.7
<50	180 (95%)	123 (95%)	57 (93%)
≥50	10 (5.3%)	6 (4.7%)	4 (6.6%)
pT				0.9
T1-T2	23 (12%)	16 (12%)	7 (11%)
T3-T4	167 (88%)	113 (88%)	54 (89%)
pN				0.015
N0	77 (41%)	60 (47%)	17 (28%)
N1-N2	113 (59%)	69 (53%)	44 (72%)
ACT				0.9
Adjuvant Chemotherapy	48 (25%)	33 (26%)	15 (25%)
Observation	142 (75%)	96 (74%)	46 (75%)
Oxaliplatin.History	50 (26%)	39 (30%)	11 (18%)	0.075
Postop.Complication	23 (12%)	10 (7.8%)	13 (21%)	0.007
BRAF.V600E				>0.9
BRAF wt	188 (99%)	127 (98%)	61 (100%)
BRAF V600E	2 (1.1%)	2 (1.6%)	0 (0%)
RAS				0.5
RAS wt	106 (56%)	74 (57%)	32 (52%)
RAS mut	84 (44%)	55 (43%)	29 (48%)
MSI				0.5
MSS	188 (99%)	128 (99%)	60 (98%)
MSI-High	2 (1.1%)	1 (0.8%)	1 (1.6%)
1Median (Range); n (%)
2Wilcoxon rank sum test; Pearson's Chi-squared test; Fisher's exact test

save_as_docx(fit1, path= "~/Downloads/merged_table.docx")

#ctDNA Detection Rates by Window and Stages

#ctDNA at MRD Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
circ_data$Stage <- factor(circ_data$Stage, levels=c("IV"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.MRD == "POSITIVE", by=list(circ_data$Stage), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.MRD, by=list(circ_data$Stage), FUN=length)
combined_data <- data.frame(
  Stage = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.MRD == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Stage = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#ctDNA at Surveillance Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels=c("NEGATIVE","POSITIVE"))
circ_data <- subset(circ_data, ctDNA.Surveillance %in% c("NEGATIVE", "POSITIVE"))
circ_data$Stage <- factor(circ_data$Stage, levels=c("IV"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.Surveillance == "POSITIVE", by=list(circ_data$Stage), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.Surveillance, by=list(circ_data$Stage), FUN=length)
combined_data <- data.frame(
  Stage = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.Surveillance == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Stage = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#ctDNA Detection Rates by Window and Synchronicity

#ctDNA at MRD Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.MRD == "POSITIVE", by=list(circ_data$Mets.Type), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.MRD, by=list(circ_data$Mets.Type), FUN=length)
combined_data <- data.frame(
  Mets.Type = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.MRD == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Mets.Type = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#ctDNA at Surveillance Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels=c("NEGATIVE","POSITIVE"))
circ_data <- subset(circ_data, ctDNA.Surveillance %in% c("NEGATIVE", "POSITIVE"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.Surveillance == "POSITIVE", by=list(circ_data$Mets.Type), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.Surveillance, by=list(circ_data$Mets.Type), FUN=length)
combined_data <- data.frame(
  Mets.Type = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.Surveillance == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Mets.Type = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#Barplot with Detection Rates at the MRD and Surveillance Windows by Synchronicity

#Detection rate at the Surveillance Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
contingency_table <- table(circ_data$Mets.Type, circ_data$ctDNA.MRD)
chi_square_test <- chisq.test(contingency_table)
print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 10.707, df = 1, p-value = 0.001067

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.0009273
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.507768 5.977774
sample estimates:
odds ratio 
  2.987004 

print(contingency_table)

              
               Negative Positive
  Metachronous       96       30
  Synchronous        33       31

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "Patients with CLM - ctDNA MRD Window", 
       x = "Synchronicity", 
       y = "Patients (%)", 
       fill = "ctDNA MRD Window",
       caption = paste("Chi-squared test p-value: ", format.pval(chi_square_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("Negative" = "blue", "Positive" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Recurrence label size

#Detection rate at the Surveillance Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]

circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
contingency_table <- table(circ_data$Mets.Type, circ_data$ctDNA.Surveillance)
chi_square_test <- chisq.test(contingency_table)
print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 4.4449, df = 1, p-value = 0.03501

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.03232
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.048856 4.688132
sample estimates:
odds ratio 
  2.210114 

print(contingency_table)

              
               Negative Positive
  Metachronous       72       30
  Synchronous        27       25

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "Patients with CLM - ctDNA Surveillance Window", 
       x = "Synchronicity", 
       y = "Patients (%)", 
       fill = "ctDNA Surveillance Window",
       caption = paste("Chi-squared test p-value: ", format.pval(chi_square_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("Negative" = "blue", "Positive" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Recurrence label size

#DFS by ctDNA at the MRD Window - Landmark MRD timepoint

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$DFS.MRD.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event)~ctDNA.MRD, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) ~ 
    ctDNA.MRD, data = circ_data)

                     n events median 0.95LCL 0.95UCL
ctDNA.MRD=NEGATIVE 129     45     NA      32      NA
ctDNA.MRD=POSITIVE  60     53   2.56       2     4.6

event_summary <- circ_data %>%
  group_by(ctDNA.MRD) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA MRD window", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(24, 30, 36))

Call: survfit(formula = surv_object ~ ctDNA.MRD, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     41      42    0.645  0.0456        0.548        0.726
   30     29       2    0.613  0.0486        0.511        0.700
   36     15       1    0.592  0.0513        0.484        0.684

                ctDNA.MRD=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24      4      53    0.108  0.0409       0.0445        0.203
   30      3       0    0.108  0.0409       0.0445        0.203
   36      2       0    0.108  0.0409       0.0445        0.203

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD, data = circ_data)

  n= 189, number of events= 98 

                    coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.MRDPOSITIVE 1.7470    5.7371   0.2092 8.352   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRDPOSITIVE     5.737     0.1743     3.808     8.645

Concordance= 0.701  (se = 0.021 )
Likelihood ratio test= 65.87  on 1 df,   p=5e-16
Wald test            = 69.75  on 1 df,   p=<2e-16
Score (logrank) test = 86.49  on 1 df,   p=<2e-16

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 5.74 (3.81-8.64); p = 0"

#Fisher test for DFS percentages at 24, 30, and 36 months
dfs_times <- c(24, 30, 36)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$DFS.MRD.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$DFS.MRD.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ctDNA.MRD == "NEGATIVE")
  pos_total <- sum(circ_data$ctDNA.MRD == "POSITIVE")
  
  neg_surv <- neg_total - sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$DFS.Event == 1 & circ_data$DFS.MRD.months < time)
  pos_surv <- pos_total - sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$DFS.Event == 1 & circ_data$DFS.MRD.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months p-value at 30 months p-value at 36 months 
        2.248794e-13         8.823505e-13         1.412511e-12 

#OS by ctDNA at the MRD Window - Landmark MRD timepoint

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$OS.MRD.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event)~ctDNA.MRD, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event) ~ 
    ctDNA.MRD, data = circ_data)

                     n events median 0.95LCL 0.95UCL
ctDNA.MRD=NEGATIVE 129      5     NA      NA      NA
ctDNA.MRD=POSITIVE  61     11   43.4      NA      NA

event_summary <- circ_data %>%
  group_by(ctDNA.MRD) %>%
  summarise(
    Total = n(),
    Events = sum(OS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - ctDNA MRD window | All stages", ylab= "Overall Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(24, 30, 36))

Call: survfit(formula = surv_object ~ ctDNA.MRD, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     71       2    0.984  0.0111        0.938        0.996
   30     48       0    0.984  0.0111        0.938        0.996
   36     26       2    0.943  0.0306        0.841        0.980

                ctDNA.MRD=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     23       7    0.836  0.0601        0.674        0.921
   30     12       3    0.717  0.0819        0.521        0.844
   36      6       0    0.717  0.0819        0.521        0.844

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD, data = circ_data)

  n= 190, number of events= 16 

                    coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.MRDPOSITIVE 1.8630    6.4429   0.5441 3.424 0.000618 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRDPOSITIVE     6.443     0.1552     2.218     18.72

Concordance= 0.726  (se = 0.064 )
Likelihood ratio test= 12.93  on 1 df,   p=3e-04
Wald test            = 11.72  on 1 df,   p=6e-04
Score (logrank) test = 15.32  on 1 df,   p=9e-05

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 4))
print(label_text)

[1] "HR = 6.44 (2.22-18.72); p = 6e-04"

#Fisher test for OS percentages at 24, 30, and 36 months
os_times <- c(24, 30, 36)
p_values <- sapply(os_times, function(time) {
  neg_count <- sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.MRD.months >= time & circ_data$OS.Event == 0)
  pos_count <- sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.MRD.months >= time & circ_data$OS.Event == 0)
  neg_total <- sum(circ_data$ctDNA.MRD == "NEGATIVE")
  pos_total <- sum(circ_data$ctDNA.MRD == "POSITIVE")
  
  neg_surv <- neg_total - sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.Event == 1 & circ_data$OS.MRD.months < time)
  pos_surv <- pos_total - sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.Event == 1 & circ_data$OS.MRD.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", os_times, " months")
print(p_values)

p-value at 24 months p-value at 30 months p-value at 36 months 
        0.0054510839         0.0002468767         0.0020570349 

#DFS by ctDNA at the MRD Window & Synchronicity

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.MRD == "NEGATIVE" & Mets.Type == "Synchronous" ~ 1,
    ctDNA.MRD == "NEGATIVE" & Mets.Type == "Metachronous" ~ 2,
    ctDNA.MRD == "POSITIVE" & Mets.Type == "Synchronous" ~ 3,
    ctDNA.MRD == "POSITIVE" & Mets.Type == "Metachronous" ~ 4
  ))

survfit(Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event)~ctDNA.Dynamics, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) ~ 
    ctDNA.Dynamics, data = circ_data)

                  n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1 33      8     NA      NA      NA
ctDNA.Dynamics=2 96     37     NA   26.91      NA
ctDNA.Dynamics=3 31     27   4.14    1.35    7.36
ctDNA.Dynamics=4 30     27   2.14    1.77    5.26

event_summary <- circ_data %>%
  group_by(ctDNA.Dynamics) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","green","purple","red"), title="DFS - ctDNA MRD window & Synchronicity", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative & Synchronous", "ctDNA Negative & Metachronous", "ctDNA Positive & Synchronous", "ctDNA Positive & Metachronous"), legend.title="")

summary(KM_curve, times= c(24, 30, 36))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Dynamics=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     15       7    0.772  0.0766        0.578        0.885
   30     10       1    0.721  0.0871        0.509        0.853
   36      6       0    0.721  0.0871        0.509        0.853

                ctDNA.Dynamics=2 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     26      35    0.600  0.0551        0.484        0.698
   30     19       1    0.576  0.0579        0.455        0.680
   36      9       1    0.546  0.0623        0.416        0.658

                ctDNA.Dynamics=3 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24      2      27    0.129  0.0602       0.0407         0.27
   30      2       0    0.129  0.0602       0.0407         0.27
   36      1       0    0.129  0.0602       0.0407         0.27

                ctDNA.Dynamics=4 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24      2      27     0.08  0.0527       0.0149         0.22
   30      1       0     0.08  0.0527       0.0149         0.22
   36      1       0     0.08  0.0527       0.0149         0.22

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1", "2", "3", "4"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

  n= 190, number of events= 99 

                  coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.Dynamics2 0.5518    1.7364   0.3901 1.415    0.157    
ctDNA.Dynamics3 2.0917    8.0983   0.4056 5.157 2.51e-07 ***
ctDNA.Dynamics4 2.3011    9.9847   0.4078 5.643 1.67e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                exp(coef) exp(-coef) lower .95 upper .95
ctDNA.Dynamics2     1.736     0.5759    0.8083      3.73
ctDNA.Dynamics3     8.098     0.1235    3.6572     17.93
ctDNA.Dynamics4     9.985     0.1002    4.4899     22.20

Concordance= 0.722  (se = 0.022 )
Likelihood ratio test= 70.35  on 3 df,   p=4e-15
Wald test            = 72.18  on 3 df,   p=1e-15
Score (logrank) test = 91.11  on 3 df,   p=<2e-16

cox_fit_summary <- summary(cox_fit)

#Repeat analysis to compare Metachronous pts by ctDNA
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.MRD == "NEGATIVE" & Mets.Type == "Synchronous" ~ 1,
    ctDNA.MRD == "NEGATIVE" & Mets.Type == "Metachronous" ~ 2,
    ctDNA.MRD == "POSITIVE" & Mets.Type == "Synchronous" ~ 3,
    ctDNA.MRD == "POSITIVE" & Mets.Type == "Metachronous" ~ 4
  ))

survfit(Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event)~ctDNA.Dynamics, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) ~ 
    ctDNA.Dynamics, data = circ_data)

                  n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1 33      8     NA      NA      NA
ctDNA.Dynamics=2 96     37     NA   26.91      NA
ctDNA.Dynamics=3 31     27   4.14    1.35    7.36
ctDNA.Dynamics=4 30     27   2.14    1.77    5.26

event_summary <- circ_data %>%
  group_by(ctDNA.Dynamics) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","green","purple","red"), title="DFS - ctDNA MRD window & Synchronicity", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative & Synchronous", "ctDNA Negative & Metachronous", "ctDNA Positive & Synchronous", "ctDNA Positive & Metachronous"), legend.title="")

summary(KM_curve, times= c(24, 30, 36))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Dynamics=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     15       7    0.772  0.0766        0.578        0.885
   30     10       1    0.721  0.0871        0.509        0.853
   36      6       0    0.721  0.0871        0.509        0.853

                ctDNA.Dynamics=2 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     26      35    0.600  0.0551        0.484        0.698
   30     19       1    0.576  0.0579        0.455        0.680
   36      9       1    0.546  0.0623        0.416        0.658

                ctDNA.Dynamics=3 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24      2      27    0.129  0.0602       0.0407         0.27
   30      2       0    0.129  0.0602       0.0407         0.27
   36      1       0    0.129  0.0602       0.0407         0.27

                ctDNA.Dynamics=4 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24      2      27     0.08  0.0527       0.0149         0.22
   30      1       0     0.08  0.0527       0.0149         0.22
   36      1       0     0.08  0.0527       0.0149         0.22

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("3", "4", "1", "2"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

  n= 190, number of events= 99 

                   coef exp(coef) se(coef)      z Pr(>|z|)    
ctDNA.Dynamics4  0.2094    1.2329   0.2731  0.767    0.443    
ctDNA.Dynamics1 -2.0917    0.1235   0.4056 -5.157 2.51e-07 ***
ctDNA.Dynamics2 -1.5398    0.2144   0.2568 -5.995 2.03e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                exp(coef) exp(-coef) lower .95 upper .95
ctDNA.Dynamics4    1.2329     0.8111   0.72196    2.1056
ctDNA.Dynamics1    0.1235     8.0983   0.05576    0.2734
ctDNA.Dynamics2    0.2144     4.6638   0.12961    0.3547

Concordance= 0.722  (se = 0.022 )
Likelihood ratio test= 70.35  on 3 df,   p=4e-15
Wald test            = 72.18  on 3 df,   p=1e-15
Score (logrank) test = 91.11  on 3 df,   p=<2e-16

cox_fit_summary <- summary(cox_fit)

# Function to perform Fisher exact test for specified groups and time points
fisher_test <- function(group1, group2, time_points) {
  p_values <- sapply(time_points, function(time) {
    group1_count <- sum(circ_data$ctDNA.Dynamics == group1 & circ_data$DFS.MRD.months >= time & circ_data$DFS.Event == 0)
    group2_count <- sum(circ_data$ctDNA.Dynamics == group2 & circ_data$DFS.MRD.months >= time & circ_data$DFS.Event == 0)
    group1_total <- sum(circ_data$ctDNA.Dynamics == group1)
    group2_total <- sum(circ_data$ctDNA.Dynamics == group2)
    
    group1_surv <- group1_total - sum(circ_data$ctDNA.Dynamics == group1 & circ_data$DFS.Event == 1 & circ_data$DFS.MRD.months < time)
    group2_surv <- group2_total - sum(circ_data$ctDNA.Dynamics == group2 & circ_data$DFS.Event == 1 & circ_data$DFS.MRD.months < time)
    
    surv_matrix <- matrix(c(group1_surv, group2_surv, group1_total - group1_surv, group2_total - group2_surv), nrow = 2)
    test_result <- fisher.test(surv_matrix)
    return(test_result$p.value)
  })
  names(p_values) <- paste0("p-value at ", time_points, " months")
  return(p_values)
}

dfs_times <- c(24, 30, 36)
p_values_1_vs_2 <- fisher_test("1", "2", dfs_times)
p_values_3_vs_4 <- fisher_test("3", "4", dfs_times)
# Print p-values for each comparison
print("P-values for ctDNA.Dynamics==1 vs ctDNA.Dynamics==2")

[1] "P-values for ctDNA.Dynamics==1 vs ctDNA.Dynamics==2"

print(p_values_1_vs_2)

p-value at 24 months p-value at 30 months p-value at 36 months 
           0.1334454            0.2041931            0.2034219 

print("P-values for ctDNA.Dynamics==3 vs ctDNA.Dynamics==4")

[1] "P-values for ctDNA.Dynamics==3 vs ctDNA.Dynamics==4"

print(p_values_3_vs_4)

p-value at 24 months p-value at 30 months p-value at 36 months 
                   1                    1                    1 

#Multivariate cox regression at MRD Window for DFS - Landmark MRD timepoint

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$DFS.MRD.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"), labels = c("Right-sided", "Left-sided"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"), labels = c("Observation", "Chemotherapy"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ctDNA.MRD + Age.Group + PrimSite + NumLiverMetsGroup + SizeLiverMetsmmGroup + pT + pN + ACT + RAS + Mets.Type + Postop.Complication + Oxaliplatin.History, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for DFS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#Multivariate cox regression at MRD Window for DFS - ctDNA MRD Positive Landmark MRD timepoint

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$DFS.MRD.months>=0,]
circ_data <- circ_data[circ_data$ctDNA.MRD=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"), labels = c("Right-sided", "Left-sided"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"), labels = c("Observation", "Chemotherapy"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ Age.Group + PrimSite + NumLiverMetsGroup + SizeLiverMetsmmGroup + pT + pN + ACT + RAS + Mets.Type + Postop.Complication + Oxaliplatin.History, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for DFS - ctDNA MRD Positive", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#Multivariate cox regression at MRD Window for DFS - ctDNA MRD Negative Landmark MRD timepoint

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$DFS.MRD.months>=0,]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"), labels = c("Right-sided", "Left-sided"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"), labels = c("Observation", "Chemotherapy"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ Age.Group + PrimSite + NumLiverMetsGroup + SizeLiverMetsmmGroup + pT + pN + ACT + RAS + Mets.Type + Postop.Complication + Oxaliplatin.History, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for DFS - ctDNA MRD Negative", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#DFS by ACT treatment in MRD negative

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 92     33     NA    25.6      NA
ACT=TRUE  33      8     NA      NA      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000      25.0000      31.0000       0.6222       0.0566       0.5011       0.7218 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000      15.0000       8.0000       0.7231       0.0859       0.5143       0.8538 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 125, number of events= 41 

           coef exp(coef) se(coef)     z Pr(>|z|)
ACTFALSE 0.5575    1.7463   0.3948 1.412    0.158

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE     1.746     0.5727    0.8054     3.786

Concordance= 0.553  (se = 0.032 )
Likelihood ratio test= 2.22  on 1 df,   p=0.1
Wald test            = 1.99  on 1 df,   p=0.2
Score (logrank) test = 2.05  on 1 df,   p=0.2

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 1.75 (0.81-3.79); p = 0.158"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
           0.3844277 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + SizeLiverMetsmmGroup + RAS + Mets.Type + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + 
    SizeLiverMetsmmGroup + RAS + Mets.Type + Oxaliplatin.History, 
    data = circ_data)

  n= 125, number of events= 41 

                           coef exp(coef) se(coef)      z Pr(>|z|)
ACTTRUE                 -0.3844    0.6809   0.4300 -0.894    0.371
GenderMale               0.3311    1.3925   0.3938  0.841    0.400
Age.Group>70             0.2674    1.3066   0.3279  0.815    0.415
NumLiverMetsGroup≥2      0.2951    1.3432   0.3658  0.807    0.420
SizeLiverMetsmmGroup≥50  0.4474    1.5642   0.6244  0.716    0.474
RASMutant                0.3107    1.3644   0.3582  0.867    0.386
Mets.TypeSynchronous    -0.1946    0.8232   0.4471 -0.435    0.663
Oxaliplatin.HistoryYes   0.2572    1.2934   0.3778  0.681    0.496

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                    0.6809     1.4687    0.2931     1.581
GenderMale                 1.3925     0.7181    0.6436     3.013
Age.Group>70               1.3066     0.7654    0.6871     2.485
NumLiverMetsGroup≥2        1.3432     0.7445    0.6559     2.751
SizeLiverMetsmmGroup≥50    1.5642     0.6393    0.4600     5.319
RASMutant                  1.3644     0.7329    0.6761     2.753
Mets.TypeSynchronous       0.8232     1.2148    0.3427     1.977
Oxaliplatin.HistoryYes     1.2934     0.7732    0.6168     2.712

Concordance= 0.638  (se = 0.042 )
Likelihood ratio test= 6.52  on 8 df,   p=0.6
Wald test            = 6.23  on 8 df,   p=0.6
Score (logrank) test = 6.37  on 8 df,   p=0.6

#DFS by ACT treatment in MRD positive

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="POSITIVE",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 31     29   1.45    0.76    3.32
ACT=TRUE  15     10  12.92    4.11      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Positive ACT vs Observation", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(12, 24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     12.0000       3.0000      27.0000       0.1290       0.0602       0.0407       0.2698 

                ACT=TRUE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12      9       6    0.600   0.126        0.318        0.797
   24      4       4    0.333   0.122        0.122        0.564

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 46, number of events= 39 

           coef exp(coef) se(coef)     z Pr(>|z|)    
ACTFALSE 1.3838    3.9899   0.3842 3.602 0.000316 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE      3.99     0.2506     1.879     8.472

Concordance= 0.675  (se = 0.03 )
Likelihood ratio test= 15  on 1 df,   p=1e-04
Wald test            = 12.97  on 1 df,   p=3e-04
Score (logrank) test = 14.61  on 1 df,   p=1e-04

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 3.99 (1.88-8.47); p = 0"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
          0.02910779 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="POSITIVE",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + SizeLiverMetsmmGroup + RAS + Mets.Type + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + 
    SizeLiverMetsmmGroup + RAS + Mets.Type + Oxaliplatin.History, 
    data = circ_data)

  n= 46, number of events= 39 

                            coef exp(coef) se(coef)      z Pr(>|z|)    
ACTTRUE                 -2.68676   0.06810  0.67360 -3.989 6.65e-05 ***
GenderMale              -1.04107   0.35308  0.51890 -2.006  0.04483 *  
Age.Group>70            -0.01813   0.98203  0.39282 -0.046  0.96318    
NumLiverMetsGroup≥2     -0.94250   0.38965  0.43212 -2.181  0.02918 *  
SizeLiverMetsmmGroup≥50  1.96295   7.12027  0.91333  2.149  0.03162 *  
RASMutant                1.01046   2.74687  0.36487  2.769  0.00562 ** 
Mets.TypeSynchronous     0.62120   1.86116  0.44165  1.407  0.15956    
Oxaliplatin.HistoryYes   0.81261   2.25379  0.57133  1.422  0.15493    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                    0.0681    14.6840   0.01819    0.2550
GenderMale                 0.3531     2.8322   0.12770    0.9762
Age.Group>70               0.9820     1.0183   0.45473    2.1208
NumLiverMetsGroup≥2        0.3897     2.5664   0.16705    0.9089
SizeLiverMetsmmGroup≥50    7.1203     0.1404   1.18869   42.6504
RASMutant                  2.7469     0.3641   1.34357    5.6159
Mets.TypeSynchronous       1.8612     0.5373   0.78316    4.4230
Oxaliplatin.HistoryYes     2.2538     0.4437   0.73552    6.9061

Concordance= 0.785  (se = 0.033 )
Likelihood ratio test= 35.06  on 8 df,   p=3e-05
Wald test            = 28.56  on 8 df,   p=4e-04
Score (logrank) test = 35.96  on 8 df,   p=2e-05

#DFS by ACT treatment in MRD negative - Pathological Stage T3-T4

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$pT=="T3-T4",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 79     32     NA    20.3      NA
ACT=TRUE  31      8     NA      NA      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation | Pathological Stage T3-T4", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000      20.0000      30.0000       0.5730       0.0626       0.4412       0.6846 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000      13.0000       8.0000       0.7009       0.0921       0.4803       0.8418 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 110, number of events= 40 

           coef exp(coef) se(coef)    z Pr(>|z|)
ACTFALSE 0.6379    1.8925   0.3961 1.61    0.107

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE     1.892     0.5284    0.8707     4.113

Concordance= 0.563  (se = 0.033 )
Likelihood ratio test= 2.92  on 1 df,   p=0.09
Wald test            = 2.59  on 1 df,   p=0.1
Score (logrank) test = 2.68  on 1 df,   p=0.1

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 1.89 (0.87-4.11); p = 0.107"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
            0.270071 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$pT=="T3-T4",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + 
    SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data = circ_data)

  n= 110, number of events= 40 

                            coef exp(coef) se(coef)      z Pr(>|z|)  
ACTTRUE                 -0.51065   0.60011  0.41623 -1.227   0.2199  
GenderMale               0.35267   1.42286  0.38584  0.914   0.3607  
Age.Group>70             0.04548   1.04653  0.34628  0.131   0.8955  
NumLiverMetsGroup≥2      0.27376   1.31489  0.37543  0.729   0.4659  
SizeLiverMetsmmGroup≥50  1.32184   3.75032  0.67677  1.953   0.0508 .
RASMutant                0.59748   1.81753  0.36546  1.635   0.1021  
Oxaliplatin.HistoryYes   0.31654   1.37237  0.37495  0.844   0.3985  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                    0.6001     1.6664    0.2654     1.357
GenderMale                 1.4229     0.7028    0.6679     3.031
Age.Group>70               1.0465     0.9555    0.5309     2.063
NumLiverMetsGroup≥2        1.3149     0.7605    0.6300     2.744
SizeLiverMetsmmGroup≥50    3.7503     0.2666    0.9954    14.130
RASMutant                  1.8175     0.5502    0.8880     3.720
Oxaliplatin.HistoryYes     1.3724     0.7287    0.6581     2.862

Concordance= 0.651  (se = 0.042 )
Likelihood ratio test= 10.01  on 7 df,   p=0.2
Wald test            = 10.52  on 7 df,   p=0.2
Score (logrank) test = 11.3  on 7 df,   p=0.1

#DFS by ACT treatment in MRD positive - Pathological Stage T3-T4

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="POSITIVE",]
circ_data <- circ_data[circ_data$pT=="T3-T4",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 26     24   1.57    1.02    4.11
ACT=TRUE  14      9  13.74    9.17      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Positive ACT vs Observation | Pathological Stage T3-T4", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(12, 24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     12.0000       3.0000      22.0000       0.1538       0.0708       0.0483       0.3146 

                ACT=TRUE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12      9       5    0.643   0.128        0.343        0.833
   24      4       4    0.357   0.128        0.130        0.594

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 40, number of events= 33 

           coef exp(coef) se(coef)     z Pr(>|z|)    
ACTFALSE 1.3686    3.9297   0.4093 3.344 0.000826 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE      3.93     0.2545     1.762     8.765

Concordance= 0.678  (se = 0.032 )
Likelihood ratio test= 12.74  on 1 df,   p=4e-04
Wald test            = 11.18  on 1 df,   p=8e-04
Score (logrank) test = 12.57  on 1 df,   p=4e-04

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 3.93 (1.76-8.76); p = 0.001"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
          0.03927147 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="POSITIVE",]
circ_data <- circ_data[circ_data$pT=="T3-T4",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + SizeLiverMetsmmGroup + RAS + Mets.Type + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + 
    SizeLiverMetsmmGroup + RAS + Mets.Type + Oxaliplatin.History, 
    data = circ_data)

  n= 40, number of events= 33 

                            coef exp(coef) se(coef)      z Pr(>|z|)    
ACTTRUE                 -2.57134   0.07643  0.70007 -3.673  0.00024 ***
GenderMale              -0.79867   0.44993  0.57010 -1.401  0.16123    
Age.Group>70             0.07987   1.08315  0.43226  0.185  0.85340    
NumLiverMetsGroup≥2     -0.92287   0.39738  0.47471 -1.944  0.05189 .  
SizeLiverMetsmmGroup≥50  2.22629   9.26542  1.21592  1.831  0.06711 .  
RASMutant                1.03251   2.80811  0.39021  2.646  0.00814 ** 
Mets.TypeSynchronous     0.66841   1.95113  0.48829  1.369  0.17104    
Oxaliplatin.HistoryYes   0.54817   1.73008  0.67452  0.813  0.41640    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                   0.07643    13.0833   0.01938    0.3014
GenderMale                0.44993     2.2226   0.14719    1.3753
Age.Group>70              1.08315     0.9232   0.46424    2.5271
NumLiverMetsGroup≥2       0.39738     2.5165   0.15672    1.0076
SizeLiverMetsmmGroup≥50   9.26542     0.1079   0.85481  100.4298
RASMutant                 2.80811     0.3561   1.30697    6.0334
Mets.TypeSynchronous      1.95113     0.5125   0.74929    5.0806
Oxaliplatin.HistoryYes    1.73008     0.5780   0.46123    6.4896

Concordance= 0.776  (se = 0.035 )
Likelihood ratio test= 26.85  on 8 df,   p=7e-04
Wald test            = 22.15  on 8 df,   p=0.005
Score (logrank) test = 26.61  on 8 df,   p=8e-04

#DFS by ACT treatment in MRD negative - Synchronous metastatic site

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$Mets.Type=="Synchronous",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 16      5     NA    9.01      NA
ACT=TRUE  17      3     NA      NA      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation | Synchronous Met Site", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000        5.000        5.000        0.661        0.125        0.364        0.844 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000        9.000        3.000        0.776        0.116        0.447        0.924 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 33, number of events= 8 

           coef exp(coef) se(coef)     z Pr(>|z|)
ACTFALSE 0.8259    2.2838   0.7334 1.126     0.26

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE     2.284     0.4379    0.5425     9.615

Concordance= 0.634  (se = 0.082 )
Likelihood ratio test= 1.32  on 1 df,   p=0.3
Wald test            = 1.27  on 1 df,   p=0.3
Score (logrank) test = 1.34  on 1 df,   p=0.2

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 2.28 (0.54-9.62); p = 0.26"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
           0.4384333 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$Mets.Type=="Synchronous",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + 
    SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data = circ_data)

  n= 33, number of events= 8 

                           coef exp(coef) se(coef)      z Pr(>|z|)  
ACTTRUE                 -1.5633    0.2094   0.8370 -1.868   0.0618 .
GenderMale              -1.1347    0.3215   0.9708 -1.169   0.2425  
Age.Group>70            -0.6898    0.5017   0.8681 -0.795   0.4268  
NumLiverMetsGroup≥2      1.0846    2.9583   0.9082  1.194   0.2324  
SizeLiverMetsmmGroup≥50  1.7630    5.8302   1.3825  1.275   0.2022  
RASMutant                1.5951    4.9288   1.0234  1.559   0.1191  
Oxaliplatin.HistoryYes       NA        NA   0.0000     NA       NA  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                    0.2094     4.7744   0.04061     1.080
GenderMale                 0.3215     3.1102   0.04796     2.155
Age.Group>70               0.5017     1.9933   0.09151     2.750
NumLiverMetsGroup≥2        2.9583     0.3380   0.49889    17.543
SizeLiverMetsmmGroup≥50    5.8302     0.1715   0.38808    87.588
RASMutant                  4.9288     0.2029   0.66321    36.630
Oxaliplatin.HistoryYes         NA         NA        NA        NA

Concordance= 0.824  (se = 0.07 )
Likelihood ratio test= 10.15  on 6 df,   p=0.1
Wald test            = 7.56  on 6 df,   p=0.3
Score (logrank) test = 9.31  on 6 df,   p=0.2

#DFS by ACT treatment in MRD negative - Metachronous metastatic site

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$Mets.Type=="Metachronous",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 76     28     NA    25.6      NA
ACT=TRUE  16      5     NA    14.6      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation | Metachronous Met Site", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000      20.0000      26.0000       0.6142       0.0632       0.4784       0.7246 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000        6.000        5.000        0.675        0.121        0.384        0.851 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 92, number of events= 33 

           coef exp(coef) se(coef)     z Pr(>|z|)
ACTFALSE 0.2494    1.2833   0.4859 0.513    0.608

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE     1.283     0.7792    0.4951     3.326

Concordance= 0.514  (se = 0.033 )
Likelihood ratio test= 0.28  on 1 df,   p=0.6
Wald test            = 0.26  on 1 df,   p=0.6
Score (logrank) test = 0.26  on 1 df,   p=0.6

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 1.28 (0.5-3.33); p = 0.608"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
                   1 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$Mets.Type=="Metachronous",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + 
    SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data = circ_data)

  n= 92, number of events= 33 

                            coef exp(coef) se(coef)      z Pr(>|z|)  
ACTTRUE                 -0.21989   0.80261  0.50394 -0.436   0.6626  
GenderMale               0.83345   2.30124  0.47385  1.759   0.0786 .
Age.Group>70             0.49048   1.63310  0.37182  1.319   0.1871  
NumLiverMetsGroup≥2      0.05637   1.05799  0.42226  0.133   0.8938  
SizeLiverMetsmmGroup≥50  0.39198   1.47990  0.75395  0.520   0.6031  
RASMutant                0.23766   1.26827  0.39924  0.595   0.5517  
Oxaliplatin.HistoryYes   0.41547   1.51508  0.39143  1.061   0.2885  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                    0.8026     1.2459    0.2989     2.155
GenderMale                 2.3012     0.4345    0.9091     5.825
Age.Group>70               1.6331     0.6123    0.7880     3.385
NumLiverMetsGroup≥2        1.0580     0.9452    0.4624     2.421
SizeLiverMetsmmGroup≥50    1.4799     0.6757    0.3377     6.486
RASMutant                  1.2683     0.7885    0.5799     2.774
Oxaliplatin.HistoryYes     1.5151     0.6600    0.7035     3.263

Concordance= 0.623  (se = 0.046 )
Likelihood ratio test= 6.09  on 7 df,   p=0.5
Wald test            = 5.71  on 7 df,   p=0.6
Score (logrank) test = 5.8  on 7 df,   p=0.6

#DFS by ACT treatment in MRD negative - 1 Liver Met

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$NumLiverMetsGroup=="1",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 69     22     NA      31      NA
ACT=TRUE  22      6     NA      NA      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation | 1 Liver Met", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000       21.000       20.000        0.674        0.063        0.534        0.780 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000       11.000        6.000        0.693        0.107        0.434        0.851 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 91, number of events= 28 

           coef exp(coef) se(coef)     z Pr(>|z|)
ACTFALSE 0.3365    1.4000   0.4618 0.729    0.466

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE       1.4     0.7143    0.5663     3.461

Concordance= 0.536  (se = 0.038 )
Likelihood ratio test= 0.56  on 1 df,   p=0.5
Wald test            = 0.53  on 1 df,   p=0.5
Score (logrank) test = 0.54  on 1 df,   p=0.5

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 1.4 (0.57-3.46); p = 0.466"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
                   1 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$NumLiverMetsGroup=="1",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + Mets.Type + SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + Mets.Type + 
    SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data = circ_data)

  n= 91, number of events= 28 

                            coef exp(coef) se(coef)      z Pr(>|z|)  
ACTTRUE                 -0.02522   0.97509  0.52282 -0.048   0.9615  
GenderMale               0.82189   2.27480  0.49890  1.647   0.0995 .
Age.Group>70             0.05878   1.06054  0.42162  0.139   0.8891  
Mets.TypeSynchronous    -0.38867   0.67796  0.55840 -0.696   0.4864  
SizeLiverMetsmmGroup≥50  1.26344   3.53757  0.63719  1.983   0.0474 *
RASMutant                0.82093   2.27261  0.45519  1.803   0.0713 .
Oxaliplatin.HistoryYes   0.25359   1.28864  0.50661  0.501   0.6167  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                    0.9751     1.0255    0.3500     2.717
GenderMale                 2.2748     0.4396    0.8556     6.048
Age.Group>70               1.0605     0.9429    0.4641     2.423
Mets.TypeSynchronous       0.6780     1.4750    0.2269     2.025
SizeLiverMetsmmGroup≥50    3.5376     0.2827    1.0147    12.333
RASMutant                  2.2726     0.4400    0.9312     5.546
Oxaliplatin.HistoryYes     1.2886     0.7760    0.4774     3.478

Concordance= 0.66  (se = 0.053 )
Likelihood ratio test= 8.34  on 7 df,   p=0.3
Wald test            = 9.02  on 7 df,   p=0.3
Score (logrank) test = 9.93  on 7 df,   p=0.2

#DFS by ACT treatment in MRD negative - 2≥ Liver Mets

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$NumLiverMetsGroup=="≥2",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 23     11   20.3    6.15      NA
ACT=TRUE  11      2     NA      NA      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation | ≥2 Liver Met", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000        4.000       11.000        0.458        0.122        0.218        0.669 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000        4.000        2.000        0.818        0.116        0.447        0.951 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 34, number of events= 13 

           coef exp(coef) se(coef)     z Pr(>|z|)
ACTFALSE 1.1484    3.1533   0.7692 1.493    0.135

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE     3.153     0.3171    0.6983     14.24

Concordance= 0.602  (se = 0.062 )
Likelihood ratio test= 2.84  on 1 df,   p=0.09
Wald test            = 2.23  on 1 df,   p=0.1
Score (logrank) test = 2.48  on 1 df,   p=0.1

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 3.15 (0.7-14.24); p = 0.135"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
            0.139822 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$NumLiverMetsGroup=="≥2",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + Mets.Type + RAS + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + Mets.Type + 
    RAS + Oxaliplatin.History, data = circ_data)

  n= 34, number of events= 13 

                            coef exp(coef)  se(coef)      z Pr(>|z|)
ACTTRUE                -1.350297  0.259163  0.833079 -1.621    0.105
GenderMale             -0.588374  0.555229  0.792352 -0.743    0.458
Age.Group>70            0.399882  1.491649  0.597240  0.670    0.503
Mets.TypeSynchronous   -0.006411  0.993610  0.794551 -0.008    0.994
RASMutant              -0.396021  0.672993  0.629514 -0.629    0.529
Oxaliplatin.HistoryYes  0.149209  1.160916  0.690752  0.216    0.829

                       exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                   0.2592     3.8586   0.05064     1.326
GenderMale                0.5552     1.8011   0.11750     2.624
Age.Group>70              1.4916     0.6704   0.46269     4.809
Mets.TypeSynchronous      0.9936     1.0064   0.20936     4.716
RASMutant                 0.6730     1.4859   0.19596     2.311
Oxaliplatin.HistoryYes    1.1609     0.8614   0.29980     4.495

Concordance= 0.661  (se = 0.086 )
Likelihood ratio test= 4.56  on 6 df,   p=0.6
Wald test            = 3.86  on 6 df,   p=0.7
Score (logrank) test = 4.16  on 6 df,   p=0.7

#DFS by ACT treatment in MRD negative - Liver Mets size <50mm

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$SizeLiverMetsmmGroup=="<50",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 87     31     NA    25.6      NA
ACT=TRUE  32      7     NA      NA      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation | Liver Met <50mm", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000      23.0000      29.0000       0.6242       0.0586       0.4983       0.7269 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000      15.0000       7.0000       0.7457       0.0856       0.5317       0.8726 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 119, number of events= 38 

           coef exp(coef) se(coef)    z Pr(>|z|)  
ACTFALSE 0.6922    1.9981   0.4194 1.65   0.0989 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE     1.998     0.5005    0.8782     4.546

Concordance= 0.567  (se = 0.032 )
Likelihood ratio test= 3.12  on 1 df,   p=0.08
Wald test            = 2.72  on 1 df,   p=0.1
Score (logrank) test = 2.83  on 1 df,   p=0.09

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 2 (0.88-4.55); p = 0.099"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
           0.2667241 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$SizeLiverMetsmmGroup=="<50",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + Mets.Type + NumLiverMetsGroup + RAS + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + Mets.Type + 
    NumLiverMetsGroup + RAS + Oxaliplatin.History, data = circ_data)

  n= 119, number of events= 38 

                          coef exp(coef) se(coef)      z Pr(>|z|)
ACTTRUE                -0.5891    0.5548   0.4541 -1.297    0.195
GenderMale              0.3234    1.3819   0.4004  0.808    0.419
Age.Group>70            0.1836    1.2016   0.3400  0.540    0.589
Mets.TypeSynchronous   -0.1471    0.8632   0.4805 -0.306    0.760
NumLiverMetsGroup≥2     0.4991    1.6472   0.3807  1.311    0.190
RASMutant               0.3732    1.4523   0.3675  1.015    0.310
Oxaliplatin.HistoryYes  0.2718    1.3123   0.3910  0.695    0.487

                       exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                   0.5548     1.8023    0.2278     1.351
GenderMale                1.3819     0.7237    0.6304     3.029
Age.Group>70              1.2016     0.8323    0.6171     2.339
Mets.TypeSynchronous      0.8632     1.1584    0.3366     2.214
NumLiverMetsGroup≥2       1.6472     0.6071    0.7811     3.474
RASMutant                 1.4523     0.6886    0.7067     2.985
Oxaliplatin.HistoryYes    1.3123     0.7620    0.6099     2.824

Concordance= 0.674  (se = 0.042 )
Likelihood ratio test= 8.79  on 7 df,   p=0.3
Wald test            = 8.67  on 7 df,   p=0.3
Score (logrank) test = 8.86  on 7 df,   p=0.3

#DFS by ACT treatment in MRD positive - Liver Mets size <50mm

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="POSITIVE",]
circ_data <- circ_data[circ_data$SizeLiverMetsmmGroup=="<50",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 29     27   1.68    1.12    3.32
ACT=TRUE  15     10  12.92    4.11      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Positive ACT vs Observation | Liver Met <50mm", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(12, 24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     12.0000       3.0000      25.0000       0.1379       0.0640       0.0435       0.2861 

                ACT=TRUE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12      9       6    0.600   0.126        0.318        0.797
   24      4       4    0.333   0.122        0.122        0.564

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 44, number of events= 37 

           coef exp(coef) se(coef)     z Pr(>|z|)    
ACTFALSE 1.3506    3.8599   0.3872 3.488 0.000487 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE      3.86     0.2591     1.807     8.245

Concordance= 0.676  (se = 0.03 )
Likelihood ratio test= 13.88  on 1 df,   p=2e-04
Wald test            = 12.17  on 1 df,   p=5e-04
Score (logrank) test = 13.65  on 1 df,   p=2e-04

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 3.86 (1.81-8.24); p = 0"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
          0.03577186 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="POSITIVE",]
circ_data <- circ_data[circ_data$SizeLiverMetsmmGroup=="<50",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + RAS + Mets.Type + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + 
    RAS + Mets.Type + Oxaliplatin.History, data = circ_data)

  n= 44, number of events= 37 

                           coef exp(coef) se(coef)      z Pr(>|z|)    
ACTTRUE                -2.67407   0.06897  0.67440 -3.965 7.34e-05 ***
GenderMale             -1.04150   0.35293  0.52608 -1.980  0.04773 *  
Age.Group>70            0.02577   1.02610  0.39461  0.065  0.94794    
NumLiverMetsGroup≥2    -0.94064   0.39038  0.43209 -2.177  0.02948 *  
RASMutant               1.01063   2.74732  0.36586  2.762  0.00574 ** 
Mets.TypeSynchronous    0.63222   1.88178  0.44608  1.417  0.15640    
Oxaliplatin.HistoryYes  0.82643   2.28514  0.57967  1.426  0.15396    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                       exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                  0.06897    14.4989   0.01839    0.2587
GenderMale               0.35293     2.8335   0.12586    0.9897
Age.Group>70             1.02610     0.9746   0.47348    2.2237
NumLiverMetsGroup≥2      0.39038     2.5616   0.16738    0.9105
RASMutant                2.74732     0.3640   1.34118    5.6277
Mets.TypeSynchronous     1.88178     0.5314   0.78499    4.5110
Oxaliplatin.HistoryYes   2.28514     0.4376   0.73366    7.1176

Concordance= 0.777  (se = 0.033 )
Likelihood ratio test= 30.27  on 7 df,   p=8e-05
Wald test            = 24.34  on 7 df,   p=0.001
Score (logrank) test = 28.16  on 7 df,   p=2e-04

#DFS by ctDNA Clearance ACT-treated at 3 months

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ACT==TRUE,]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.MRD == "POSITIVE" & ctDNA.3months == "NEGATIVE" ~ 1,
    ctDNA.MRD == "POSITIVE" & ctDNA.3months == "POSITIVE" ~ 2
  ))

circ_data <- circ_data[circ_data$DFS.3mo.months>=0,]
survfit(Surv(time = circ_data$DFS.3mo.months, event = circ_data$DFS.Event)~ctDNA.Dynamics, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.3mo.months, event = circ_data$DFS.Event) ~ 
    ctDNA.Dynamics, data = circ_data)

   35 observations deleted due to missingness 
                 n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1 8      4  14.00   13.80      NA
ctDNA.Dynamics=2 5      5   3.42    3.02      NA

event_summary <- circ_data %>%
  group_by(ctDNA.Dynamics) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.3mo.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA Clearance from MRD to 3 months ACT-treated", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Clearance", "No Clearance"), legend.title="")

summary(KM_curve, times= c(6, 12, 24))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

35 observations deleted due to missingness 
                ctDNA.Dynamics=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    6      7       1    0.875   0.117        0.387        0.981
   12      7       0    0.875   0.117        0.387        0.981
   24      3       3    0.500   0.177        0.152        0.775

                ctDNA.Dynamics=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
       6.000        2.000        3.000        0.400        0.219        0.052        0.753 

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

  n= 13, number of events= 9 
   (35 observations deleted due to missingness)

                             coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.DynamicsNo Clearance  2.717    15.131    1.117 2.432    0.015 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.DynamicsNo Clearance     15.13    0.06609     1.695     135.1

Concordance= 0.736  (se = 0.053 )
Likelihood ratio test= 8.61  on 1 df,   p=0.003
Wald test            = 5.92  on 1 df,   p=0.02
Score (logrank) test = 9.83  on 1 df,   p=0.002

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 15.13 (1.69-135.08); p = 0.015"

#DFS by ctDNA Clearance ACT-treated at 6 months

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ACT==TRUE,]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.MRD == "POSITIVE" & ctDNA.6months == "NEGATIVE" ~ 1,
    ctDNA.MRD == "POSITIVE" & ctDNA.6months == "POSITIVE" ~ 2
  ))

circ_data <- circ_data[circ_data$DFS.6mo.months>=0,]
survfit(Surv(time = circ_data$DFS.6mo.months, event = circ_data$DFS.Event)~ctDNA.Dynamics, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.6mo.months, event = circ_data$DFS.Event) ~ 
    ctDNA.Dynamics, data = circ_data)

   36 observations deleted due to missingness 
                 n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1 9      5 11.696   9.429      NA
ctDNA.Dynamics=2 3      3  0.789   0.493      NA

event_summary <- circ_data %>%
  group_by(ctDNA.Dynamics) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.6mo.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA Clearance from MRD to 6 months ACT-treated", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Clearance", "No Clearance"), legend.title="")

summary(KM_curve, times= c(6, 24))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

36 observations deleted due to missingness 
                ctDNA.Dynamics=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    6      7       2    0.778   0.139        0.365        0.939
   24      2       3    0.444   0.166        0.136        0.719

                ctDNA.Dynamics=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     6.00000      1.00000      2.00000      0.33333      0.27217      0.00896      0.77415 

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

  n= 12, number of events= 8 
   (36 observations deleted due to missingness)

                             coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.DynamicsNo Clearance 1.9221    6.8355   0.9345 2.057   0.0397 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.DynamicsNo Clearance     6.836     0.1463     1.095     42.68

Concordance= 0.658  (se = 0.073 )
Likelihood ratio test= 4.19  on 1 df,   p=0.04
Wald test            = 4.23  on 1 df,   p=0.04
Score (logrank) test = 5.53  on 1 df,   p=0.02

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 6.84 (1.09-42.68); p = 0.04"

#DFS by ctDNA at the Surveillance Window - Landmark 10 weeks

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Surveillance!="",]
circ_data$DFS.months=circ_data$DFS.months-2.5
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ctDNA.Surveillance, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ctDNA.Surveillance, data = circ_data)

                             n events median 0.95LCL 0.95UCL
ctDNA.Surveillance=NEGATIVE 99     18     NA      NA      NA
ctDNA.Surveillance=POSITIVE 53     44   5.81    3.05    12.1

event_summary <- circ_data %>%
  group_by(ctDNA.Surveillance) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Surveillance, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA Surveillance window", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(24, 30, 36))

Call: survfit(formula = surv_object ~ ctDNA.Surveillance, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Surveillance=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     41      16    0.815  0.0434        0.711        0.884
   30     31       1    0.795  0.0468        0.684        0.870
   36     15       1    0.768  0.0522        0.646        0.853

                ctDNA.Surveillance=POSITIVE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000       2.0000      44.0000       0.0859       0.0501       0.0200       0.2142 

circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.Surveillance, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Surveillance, data = circ_data)

  n= 152, number of events= 62 

                              coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.SurveillancePOSITIVE  2.3644   10.6382   0.2982 7.929 2.21e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.SurveillancePOSITIVE     10.64      0.094      5.93     19.09

Concordance= 0.757  (se = 0.025 )
Likelihood ratio test= 73.09  on 1 df,   p=<2e-16
Wald test            = 62.87  on 1 df,   p=2e-15
Score (logrank) test = 91.54  on 1 df,   p=<2e-16

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 10.64 (5.93-19.09); p = 0"

#Fisher test for DFS percentages at 24, 30, and 36 months
dfs_times <- c(24, 30, 36)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ctDNA.Surveillance == "NEGATIVE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ctDNA.Surveillance == "POSITIVE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ctDNA.Surveillance == "NEGATIVE")
  pos_total <- sum(circ_data$ctDNA.Surveillance == "POSITIVE")
  
  neg_surv <- neg_total - sum(circ_data$ctDNA.Surveillance == "NEGATIVE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ctDNA.Surveillance == "POSITIVE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months p-value at 30 months p-value at 36 months 
        4.109654e-16         1.518274e-15         4.230750e-15 

#OS by ctDNA at the Surveillance Window - Landmark 10 weeks

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Surveillance!="",]
circ_data$OS.months=circ_data$OS.months-2.5
circ_data <- circ_data[circ_data$OS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$OS.months, event = circ_data$OS.Event)~ctDNA.Surveillance, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$OS.months, event = circ_data$OS.Event) ~ 
    ctDNA.Surveillance, data = circ_data)

                             n events median 0.95LCL 0.95UCL
ctDNA.Surveillance=NEGATIVE 99      3     NA      NA      NA
ctDNA.Surveillance=POSITIVE 55      6   41.8    37.3      NA

event_summary <- circ_data %>%
  group_by(ctDNA.Surveillance) %>%
  summarise(
    Total = n(),
    Events = sum(OS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$OS.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Surveillance, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - ctDNA Surveillance window", ylab= "Overall Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(24, 30, 36))

Call: survfit(formula = surv_object ~ ctDNA.Surveillance, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Surveillance=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     55       1    0.990  0.0104        0.928        0.999
   30     39       1    0.965  0.0259        0.855        0.992
   36     18       1    0.941  0.0351        0.817        0.982

                ctDNA.Surveillance=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     22       2    0.955  0.0319        0.827        0.989
   30     11       2    0.858  0.0715        0.641        0.948
   36      5       0    0.858  0.0715        0.641        0.948

circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.Surveillance, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Surveillance, data = circ_data)

  n= 154, number of events= 9 

                             coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.SurveillancePOSITIVE 1.6394    5.1523   0.7188 2.281   0.0226 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.SurveillancePOSITIVE     5.152     0.1941     1.259     21.08

Concordance= 0.678  (se = 0.094 )
Likelihood ratio test= 5.59  on 1 df,   p=0.02
Wald test            = 5.2  on 1 df,   p=0.02
Score (logrank) test = 6.33  on 1 df,   p=0.01

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 5.15 (1.26-21.08); p = 0.023"

#Fisher test for OS percentages at 24, 30, and 36 months
os_times <- c(24, 30, 36)
p_values <- sapply(os_times, function(time) {
  neg_count <- sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.months >= time & circ_data$OS.Event == 0)
  pos_count <- sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.months >= time & circ_data$OS.Event == 0)
  neg_total <- sum(circ_data$ctDNA.MRD == "NEGATIVE")
  pos_total <- sum(circ_data$ctDNA.MRD == "POSITIVE")
  
  neg_surv <- neg_total - sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.Event == 1 & circ_data$OS.months < time)
  pos_surv <- pos_total - sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.Event == 1 & circ_data$OS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", os_times, " months")
print(p_values)

p-value at 24 months p-value at 30 months p-value at 36 months 
          0.12280702           0.02179403           0.04276279 

#Multivariate cox regression at Surveillance Window for DFS - Landmark 10 weeks

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Surveillance!="",]
circ_data$DFS.months=circ_data$DFS.months-2.5
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels=c("NEGATIVE","POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"), labels = c("Right-sided", "Left-sided"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"), labels = c("Observation", "Chemotherapy"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ctDNA.Surveillance + Age.Group + PrimSite + NumLiverMetsGroup + SizeLiverMetsmmGroup + pT + pN + ACT + RAS + Mets.Type + Postop.Complication + Oxaliplatin.History, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for DFS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#OS by ctDNA at the MRD Window - pts with Radiological Recurrence

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$RFS.Event=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$OS.MRD.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event)~ctDNA.MRD, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event) ~ 
    ctDNA.MRD, data = circ_data)

                    n events median 0.95LCL 0.95UCL
ctDNA.MRD=NEGATIVE 42      2     NA      NA      NA
ctDNA.MRD=POSITIVE 53     10   43.4      NA      NA

event_summary <- circ_data %>%
  group_by(ctDNA.MRD) %>%
  summarise(
    Total = n(),
    Events = sum(OS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - Radiological Recurrence | ctDNA MRD window", ylab= "Overall Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(24, 36))

Call: survfit(formula = surv_object ~ ctDNA.MRD, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     25       0    1.000  0.0000           NA           NA
   36      8       1    0.933  0.0644        0.613         0.99

                ctDNA.MRD=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     18       6     0.82  0.0709        0.627        0.919
   36      4       3     0.68  0.0942        0.459        0.827

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD, data = circ_data)

  n= 95, number of events= 12 

                    coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.MRDPOSITIVE 1.8047    6.0781   0.7937 2.274    0.023 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRDPOSITIVE     6.078     0.1645     1.283      28.8

Concordance= 0.735  (se = 0.031 )
Likelihood ratio test= 6.88  on 1 df,   p=0.009
Wald test            = 5.17  on 1 df,   p=0.02
Score (logrank) test = 6.56  on 1 df,   p=0.01

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 6.08 (1.28-28.8); p = 0.023"

#Fisher test for OS percentages at 24 and 36 months
os_times <- c(24, 36)
p_values <- sapply(os_times, function(time) {
  neg_count <- sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.MRD.months >= time & circ_data$OS.Event == 0)
  pos_count <- sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.MRD.months >= time & circ_data$OS.Event == 0)
  neg_total <- sum(circ_data$ctDNA.MRD == "NEGATIVE")
  pos_total <- sum(circ_data$ctDNA.MRD == "POSITIVE")
  
  neg_surv <- neg_total - sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.Event == 1 & circ_data$OS.MRD.months < time)
  pos_surv <- pos_total - sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.Event == 1 & circ_data$OS.MRD.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", os_times, " months")
print(p_values)

p-value at 24 months p-value at 36 months 
          0.03245083           0.03892912 

#PRS by ctDNA at the MRD Window - pts with Radiological Recurrence

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$RFS.Event=="TRUE",]
circ_data <- circ_data[circ_data$OS.MRD.months>=0,]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]

survfit(Surv(time = circ_data$PRS.months, event = circ_data$OS.Event)~ctDNA.MRD, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$PRS.months, event = circ_data$OS.Event) ~ 
    ctDNA.MRD, data = circ_data)

                    n events median 0.95LCL 0.95UCL
ctDNA.MRD=NEGATIVE 42      2     NA    36.3      NA
ctDNA.MRD=POSITIVE 53     10   38.2      NA      NA

event_summary <- circ_data %>%
  group_by(ctDNA.MRD) %>%
  summarise(
    Total = n(),
    Events = sum(OS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$PRS.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="PRS - Radiological Recurrence | ctDNA MRD window", ylab= "Post-Recurrence Survival", xlab="Time from Radiological Recurrence (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(24, 30))

Call: survfit(formula = surv_object ~ ctDNA.MRD, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     14       1    0.938  0.0605        0.632        0.991
   30      7       0    0.938  0.0605        0.632        0.991

                ctDNA.MRD=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     13       8    0.731  0.0894        0.510        0.864
   30      6       1    0.675  0.0986        0.443        0.827

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD, data = circ_data)

  n= 95, number of events= 12 

                    coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.MRDPOSITIVE 1.4643    4.3246   0.7778 1.883   0.0597 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRDPOSITIVE     4.325     0.2312    0.9417     19.86

Concordance= 0.683  (se = 0.042 )
Likelihood ratio test= 4.67  on 1 df,   p=0.03
Wald test            = 3.54  on 1 df,   p=0.06
Score (logrank) test = 4.2  on 1 df,   p=0.04

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 4.32 (0.94-19.86); p = 0.06"

#Fisher test for PRS percentages at 24 and 30 months
prs_times <- c(24, 30)
p_values <- sapply(prs_times, function(time) {
  neg_count <- sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$PRS.months >= time & circ_data$OS.Event == 0)
  pos_count <- sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$PRS.months >= time & circ_data$OS.Event == 0)
  neg_total <- sum(circ_data$ctDNA.MRD == "NEGATIVE")
  pos_total <- sum(circ_data$ctDNA.MRD == "POSITIVE")
  
  neg_surv <- neg_total - sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.Event == 1 & circ_data$PRS.months < time)
  pos_surv <- pos_total - sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.Event == 1 & circ_data$PRS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", prs_times, " months")
print(p_values)

p-value at 24 months p-value at 30 months 
          0.07279604           0.03892912 

#Percentage of ctDNA MRD Window positivity in pts with liver metastasis

rm(list = ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data %>%
  filter(LiverMets == "TRUE" & RFS.Event == "TRUE" & ctDNA.MRD != "")
circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))

positive_rate <- sum(circ_data$ctDNA.MRD == "Positive" & grepl("liver", circ_data$RelSite, ignore.case = TRUE)) / sum(circ_data$ctDNA.MRD == "Positive") * 100
positive_ci <- binconf(sum(circ_data$ctDNA.MRD == "Positive" & grepl("liver", circ_data$RelSite, ignore.case = TRUE)),
                       sum(circ_data$ctDNA.MRD == "Positive"),
                       alpha = 0.05)[c(2, 3)] * 100
negative_rate <- sum(circ_data$ctDNA.MRD == "Negative" & grepl("liver", circ_data$RelSite, ignore.case = TRUE)) / sum(circ_data$ctDNA.MRD == "Negative") * 100
negative_ci <- binconf(sum(circ_data$ctDNA.MRD == "Negative" & grepl("liver", circ_data$RelSite, ignore.case = TRUE)),
                       sum(circ_data$ctDNA.MRD == "Negative"),
                       alpha = 0.05)[c(2, 3)] * 100
data <- data.frame(
  ctDNA.MRD = c("Positive", "Negative"),
  percentage = c(positive_rate, negative_rate),
  lower_ci = c(positive_ci[1], negative_ci[1]),
  upper_ci = c(positive_ci[2], negative_ci[2])
)
cross_tab <- table(circ_data$ctDNA.MRD, grepl("liver", circ_data$RelSite, ignore.case = TRUE))
chi_test <- chisq.test(cross_tab)
p_value <- format.pval(chi_test$p.value, digits = 4)
print(data)

print(cross_tab)

          
           FALSE TRUE
  Negative    28   14
  Positive    14   39

print(chi_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  cross_tab
X-squared = 13.803, df = 1, p-value = 0.000203

barplot <- ggplot(data, aes(x = ctDNA.MRD, y = percentage, fill = ctDNA.MRD)) +
  geom_bar(stat = "identity") +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2) +
  geom_text(aes(label = paste0(round(percentage, 1), "%")), vjust = -0.5) +
  labs(
    x = "ctDNA status at the MRD status",
    y = "Proportion of patients with recurrence in the liver",
    caption = paste("Chi-squared test p-value: ", p_value)
  ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 100)) +
  scale_fill_manual(values = c("Negative" = "blue", "Positive" = "red")) +
  theme_minimal()
print(barplot)

#Percentage of ctDNA MRD Window positivity in pts with lung metastasis

rm(list = ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data %>%
  filter(LiverMets == "TRUE" & RFS.Event == "TRUE" & ctDNA.MRD != "")
circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))

positive_rate <- sum(circ_data$ctDNA.MRD == "Positive" & grepl("lung", circ_data$RelSite, ignore.case = TRUE)) / sum(circ_data$ctDNA.MRD == "Positive") * 100
positive_ci <- binconf(sum(circ_data$ctDNA.MRD == "Positive" & grepl("lung", circ_data$RelSite, ignore.case = TRUE)),
                       sum(circ_data$ctDNA.MRD == "Positive"),
                       alpha = 0.05)[c(2, 3)] * 100
negative_rate <- sum(circ_data$ctDNA.MRD == "Negative" & grepl("lung", circ_data$RelSite, ignore.case = TRUE)) / sum(circ_data$ctDNA.MRD == "Negative") * 100
negative_ci <- binconf(sum(circ_data$ctDNA.MRD == "Negative" & grepl("lung", circ_data$RelSite, ignore.case = TRUE)),
                       sum(circ_data$ctDNA.MRD == "Negative"),
                       alpha = 0.05)[c(2, 3)] * 100
data <- data.frame(
  ctDNA.MRD = c("Positive", "Negative"),
  percentage = c(positive_rate, negative_rate),
  lower_ci = c(positive_ci[1], negative_ci[1]),
  upper_ci = c(positive_ci[2], negative_ci[2])
)
cross_tab <- table(circ_data$ctDNA.MRD, grepl("lung", circ_data$RelSite, ignore.case = TRUE))
chi_test <- chisq.test(cross_tab)
p_value <- format.pval(chi_test$p.value, digits = 4)
print(data)

print(cross_tab)

          
           FALSE TRUE
  Negative    16   26
  Positive    39   14

print(chi_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  cross_tab
X-squared = 10.695, df = 1, p-value = 0.001074

barplot <- ggplot(data, aes(x = ctDNA.MRD, y = percentage, fill = ctDNA.MRD)) +
  geom_bar(stat = "identity") +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2) +
  geom_text(aes(label = paste0(round(percentage, 1), "%")), vjust = -0.5) +
  labs(
    x = "ctDNA status at the MRD status",
    y = "Proportion of patients with recurrence in the lung",
    caption = paste("Chi-squared test p-value: ", p_value)
  ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 100)) +
  scale_fill_manual(values = c("Negative" = "blue", "Positive" = "red")) +
  theme_minimal()
print(barplot)

#Barplot with Recurrence Sites (Liver vs Others) by ctDNA at the MRD Window

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$RFS.Event=="TRUE",]

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Liver <- factor(circ_data$Liver, levels = c("FALSE", "TRUE"), labels = c("Others", "Liver"))
contingency_table <- table(circ_data$ctDNA.MRD, circ_data$Liver)
chi_square_test <- chisq.test(contingency_table)
print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 18.352, df = 1, p-value = 1.836e-05

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 9.391e-06
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  2.744866 20.863460
sample estimates:
odds ratio 
  7.311049 

print(contingency_table)

          
           Others Liver
  Negative     30    12
  Positive     13    39

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "Patients with Radiological Recurrence", 
       x = "ctDNA at the MRD Window", 
       y = "Patients (%)", 
       fill = "Recurrence Site",
       caption = paste("Chi-squared test p-value: ", format.pval(chi_square_test$p.value))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("Others" = "blue", "Liver" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Recurrence label size

#Table with recurrence sites by ctDNA at the MRD Window

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$RFS.Event=="TRUE",]

circ_data_subset1 <- circ_data %>%
  select(
    Rec.Site) %>%
  mutate(
    Rec.Site = factor(Rec.Site, levels = c("Local/LN", "Liver", "Lung", "Peritoneum", "Peritoneum & Others", "Brain")))

circ_data1 <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data1 <- circ_data1[circ_data1$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$RFS.Event=="TRUE",]

circ_data_subset2 <- circ_data %>%
  select(
    Rec.Site,
    ctDNA.MRD) %>%
  mutate(
    Rec.Site = factor(Rec.Site, levels = c("Local/LN", "Liver", "Lung", "Peritoneum", "Peritoneum & Others", "Brain")),
    ctDNA.MRD = factor(ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE")))
Overall <- circ_data_subset1 %>%
  tbl_summary(
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)")) %>%
  bold_labels()
Overall

Characteristic	N = 951
Rec.Site
Local/LN	4 (4.3%)
Liver	48 (51%)
Lung	30 (32%)
Peritoneum	4 (4.3%)
Peritoneum & Others	7 (7.4%)
Brain	1 (1.1%)
Unknown	1
1 n (%)

ByctDNA_MRD <- circ_data_subset2 %>%
  tbl_summary(
    by = ctDNA.MRD, # add this line to subgroup by ctDNA.MRD
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)")) %>%
  add_p() %>%
  bold_labels()
ByctDNA_MRD

Characteristic	NEGATIVE, N = 421	POSITIVE, N = 531	p-value2
Rec.Site			<0.001
Local/LN	1 (2.4%)	3 (5.8%)
Liver	12 (29%)	36 (69%)
Lung	23 (55%)	7 (13%)
Peritoneum	3 (7.1%)	1 (1.9%)
Peritoneum & Others	3 (7.1%)	4 (7.7%)
Brain	0 (0%)	1 (1.9%)
Unknown	0	1
1 n (%)
2 Fisher’s exact test

merged_table <- tbl_merge(tbls=list(Overall, ByctDNA_MRD))
merged_table

Characteristic	Table 1	Table 2
	N = 951	NEGATIVE, N = 421	POSITIVE, N = 531	p-value2
Rec.Site				<0.001
Local/LN	4 (4.3%)	1 (2.4%)	3 (5.8%)
Liver	48 (51%)	12 (29%)	36 (69%)
Lung	30 (32%)	23 (55%)	7 (13%)
Peritoneum	4 (4.3%)	3 (7.1%)	1 (1.9%)
Peritoneum & Others	7 (7.4%)	3 (7.1%)	4 (7.7%)
Brain	1 (1.1%)	0 (0%)	1 (1.9%)
Unknown	1	0	1
1 n (%)
2 Fisher’s exact test

fit1 <- as_flex_table(
  merged_table,
  include = everything(),
  return_calls = FALSE,
  strip_md_bold = TRUE)
fit1

	Table 1	Table 2
Characteristic	N = 951	NEGATIVE, N = 421	POSITIVE, N = 531	p-value2
Rec.Site				<0.001
Local/LN	4 (4.3%)	1 (2.4%)	3 (5.8%)
Liver	48 (51%)	12 (29%)	36 (69%)
Lung	30 (32%)	23 (55%)	7 (13%)
Peritoneum	4 (4.3%)	3 (7.1%)	1 (1.9%)
Peritoneum & Others	7 (7.4%)	3 (7.1%)	4 (7.7%)
Brain	1 (1.1%)	0 (0%)	1 (1.9%)
Unknown	1	0	1
1n (%)
2Fisher's exact test

save_as_docx(fit1, path= "~/Downloads/merged_table.docx")